Plykin-like attractor in non-autonomous coupled oscillators 
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A system of two coupled non-autonomous oscillators is considered. Dynamics of complex am- 
plitudes is governed by differential equations with periodic piecewise continuous dependence of the 
coefficients on time. The Poincare map is derived explicitly. With exclusion of the overall phase, on 
which the evolution of other variables does not depend, the Poincare map is reduced to 3D mapping. 
It possesses an attractor of Plykin type located on an invariant sphere. Computer verification of 
the cone criterion confirms the hyperbolic nature of the attractor in the 3D map. Some results of 
numerical studies of the dynamics for the coupled oscillators are presented, including the attractor 
portraits, Lyapunov exponents, and the power spectral density. 

PACS numbers: 05.45.-a, 05.40.Ca 



In mathematical theory of dynamical systems 
a class of uniformly hyperbolic strange attrac- 
tors is known. In such an attractor all orbits 
are of the same saddle type, they manifest strong 
stochastic properties and allow detailed theoret- 
ical analysis. The mathematical theory was ad- 
vanced more than 40 years ago, but till now the 
hyperbolic strange attractors are regarded rather 
as purified image of deterministic chaos than as 
realistic models of complex dynamics. In text- 
books and reviews, examples of these attractors 
are traditionally represented by abstract artifi- 
cial constructions like the Plykin attractor and 
the Smale - Williams attractor. Recently, a re- 
alistic system was suggested and implemented as 
electronic device, dynamics of which in strobo- 
scopic description is associated with attractor of 
Smale - Williams type. In the present article I 
show how the dynamics related to an attractor 
of Plykin type may be obtained in coupled non- 
autonomous self-oscillators. As systems of cou- 
pled oscillators occur in many fields in physics 
and technology, it is natural to expect that the 
suggested model may be realizable, for example, 
with electronic devices, mechanical systems, ob- 
jects of laser physics and nonlinear optics. The 
systems with hyperbolic strange attractors may 
be of special interest in applications (e.g. for 
noise generators, chaos communication etc.) due 
to the intrinsic structural stability, that means in- 
sensitivity of the chaotic motions to variations of 
parameters, characteristics of elements, technical 
fluctuations etc. 
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I. INTRODUCTION 

Mathematical theory of dynamical systems introduces 
a class of uniformly hyperbolic strange attractors [ll, [3, 
3, 4, 5, 6, 7, 8, In such an attractor all orbits are of 
saddle type, and their stable and unstable manifolds do 
not touch, but can only intersect transversally. These at- 
tractors manifest strong stochastic properties and allow 
a detailed mathematical analysis. They are structurally 
stable; that means insensitivity of the structure of the 
attractors in respect to variation of functions and pa- 
rameters in the dynamical equations. Until very recent 
times, the hyperbolic strange attractors were regarded 
rather as purified image of chaos than as objects relating 
to complex dynamics of real-world systems. (See discus- 
sion of the question in Ref. 9]; also, a mechanical system 
with hyperbolic dynamics, the so-called triple linkage, 
has been considered there.) 

In textbooks and reviews, examples of the uniformly 
hyperbolic attractors are traditionally represented by 
mathematical constructions, the Plykin attractor and the 
Smale - Williams solenoid. These examples relate to 
discrete-time systems, the iterated maps. The Smale - 
Williams attractor appears in the mapping of a toroidal 
domain into itself in the state space of dimension 3 or 
more. The Plykin attractor occurs in some special map- 
ping on a sphere with four holes, or in a bounded do- 
main on a plane with three holes (Fig. 1 a) It is 
known that a variety of topologically different Plykin-like 
attractors may be constructed in finite two-dimensional 
domains with holes. One of the modifications shown in 
Fig. 1 b is of special interest for the present study and 
will be referred to as the Plykin - Newhouse attractor 

In applied disciplines, physics and technology, people 
deal more often with systems operating in continuous 
time; they are called the flows in mathematical literature. 
The procedure of passage from mapping x„j_i = f (x„) 
to a flow system is called suspension 0, H, |j, H, S 0]- 
Such a passage is possible if the map is reversible. For 
the resulting flow system the relation x„+i = f (x„) is the 
Poincare map, which in the context on non-autonomous 
systems is called sometimes the stroboscopic map. 



2 



Recently, a system was suggested and realized experi- 
mentally, in which the Poincare map possesses an attrac- 
tor of Smale - Williams type [13, It is composed 

of two non- autonomous van der Pol oscillators, which 
become active turn by turn and transfer the excitation 
each other, in such manner that the transformation of 
the phase of oscillations on a whole cycle corresponds to 
expanding circle map. Computer verification of condi- 
tions guaranteeing the hyperbolic nature of the attractor 
was performed in Ref. [l4|. (See some developments of 
the scheme in Refs. 0,E1, [II El-) 

Till now, no explicit examples were advanced for a 
Plykin type attractor to occur in a low-dimensional phys- 
ically realizable system [s^j- In Refs. [lO] and [2l| the 
authors argue in favor of existence of the Plykin-type 
attractors in the Poincare maps for a modified Lorenz 
system and for an autonomous three-dimensional system 
modeling dynamics of neuron. On the other hand, an 
explicit example of a non-autonomous flow system with 
Plykin-Newhouse attractor in the stroboscopic map has 
been advanced in the PhD thesis of Hunt [23| . The model 
of Hunt is defined by multiple expressions, distinct for 
different domains in the state space, and contains many 
artificially introduced smoothing factors. It is really hard 
to imagine that this model could be reproduced on a base 
of some physical system. 




(b) 

FIG. 1; Illustration of action on a plane for the map suggested 
in the original paper of Plykin (a) , and for the version of the 
map with the Plykin - Newhouse attractor (b) . Each of them 
may be associated with a map defined on the sphere, say, by 
means of the stereographic projection 

In the present article I show how the dynamics asso- 
ciated with attractor of Plykin type may be obtained in 
a system of coupled non- autonomous oscillators. As be- 
lieved, it opens prospects for constructing physical and 
technical systems, e.g. electron devices with the struc- 
turally stable chaotic regimes. 

In Section |TT] a sequence of continuous transformations 
is defined on a two-dimensional sphere, and a system of 
two coupled oscillators is introduced, in which the state 



evolution corresponds in some sense to those transfor- 
mations. The equations are written down for complex 
amplitudes of the oscillations. The points on the sphere 
represent the instantaneous states defined up to the over- 
all phase factor. An explicit Poincare map is derived 
that describes evolution of the state on one period of 
variation of coefficients in the non- autonomous differen- 
tial equations. With exclusion of the overall phase, on 
which the evolution of other variables does not depend, 
the Poincare map is reduced to a three-dimensional map, 
which possesses an attractor of Plykin type on an invari- 
ant sphere. In Section IIIII results of computer verifica- 
tion of the so-called cone criterion are presented confirm- 
ing the hyperbolic nature of the attractor of the three- 
dimensional map; it means also its structural stability. 
The topological type of the attractor corresponds to the 
construction of Plykin - Newhouse. In Section IIVI some 
results of numerical studies of the dynamics of the cou- 
pled oscillator system are discussed, including portraits 
of the attractor, Lyapunov exponents, power spectral 
density. In the set of equations for complex amplitudes, 
because of presence of a neutral direction in the state 
space, which is associated with the overall phase, the at- 
tractor has to be related formally to the class of partially 
hyperbolic ones d, [2^ . 



II. REPRESENTATION OF STATES ON 
A SPHERE AND EQUATIONS DESCRIBING 
DYNAMICS OF THE MODEL 

Let us start with a system of two self-oscillators with 
compensation of losses from the common energy source. 
Let the equations for the slow amplitudes a and b read 

i^(l_|a|2_|fe|2)a, 

(1) 

where /i is a positive parameter. Let us set 

b = ^e''^/2+''^ sin(6'/2), 

(2) 

a = ^e-'v/'^+^^ cos(6'/2). 

Clearly, in sustained regime of self-oscillations, the con- 
dition p = |ap + |6p = 1 has to be valid. If we consider 
states satisfying p = 1 and identify the states differing 
only in the overall phase, we can associate them with the 
points on a unit sphere (Fig. 2). Also, on the picture the 
Cartesian coordinates are shown: 

X = p cos ip sin 0, 

y — psmipsm.9, (3) 
z — pcos9. 

Via the complex amplitudes they are expressed as 

x + iy^ 2a*b, z = |ap - \b\^ . (4) 
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We intend to modify the model ([T]) in order to obtain 
a set of equations with coefficients periodically varying 
in time, in such way that in the stroboscopic description 
and in the representation on the unit sphere, the Plykin 
type attractor will occur. 

As proved by Plykin, a uniformly hyperbolic attractor 
may exist on a sphere only in presence of at least four 
holes, the areas not visiting by trajectories belonging to 
the attractor. In our construction, the holes will corre- 
spond to neighborhoods of four points A, B, C, D, with 
coordinates {x,y,z) = (±1/V2, 0, ±l/-\/2). 

Let us consider a sequence of the following continuous 
transformations on the sphere: 

• Flow down along circles of latitude, that is 
motion of the representative points on the sphere 
away from the meridians NABS and NDCS towards 
the meridians equally distant from the arcs AB and 
CD. 

• Differential rotation around z-axis with angular 
velocity depending on z linearly, in such way that 
the points B and C do not move, while the points 
A and D exchange their location. 

• Rotation of the sphere by 90° around the y-axis. 

• Flow down along circles of latitude, like at the 
first stage. 

• Inverse differential rotation around z-axis. 

• Inverse rotation by 90° around the y-axis. 

The procedure is symmetric in the sense that the op- 
erations for the stages (I) and (IV) are identical, while 
the stages (II) and (III) differ from (V) and (VI) only by 
directions of the rotations. Intuitively, it looks reason- 
able that this sequence of transformations will generate 
a flow on the sphere accompanying with formation of fil- 
aments of fine transversal structure, presence of which is 
a characteristic feature of the Plykin type attractors. 

Let us construct equations for the complex amplitudes 
to reproduce dynamics on the above stages with corre- 
sponding motion of the points on the unit square repre- 
senting the instantaneous states of the system. Duration 
of each of the six stages is accepted to be equal to a unit 
time interval. 

On a stage of flow down along circles of latitude 

we require the angular velocity of motion on the sphere 
to be proportional to sin 2ip. The simplest appropriate 
form of the differential equations is 



-i£alm(a*^6^), b = ieblm{a*^b'^) 



(5) 



Indeed, substituting 6 = e*'^/2+^^ sin(6'/2) and a = 
^-iip/2+iip cos(6'/2), after some simple transformations we 
get (fi — ^esin^ 6'sin2(/3, ^ = 0. Physically, the terms in 
the right-hand parts of ([5]) give rise to a frequency shift of 
opposite sign for two oscillators. Magnitude of the shift 
is proportional to the amplitude of a low-frequency signal 




FIG. 2: The unit sphere with marked points A, B, C, D, 
neighborhoods of which in the further construction will cor- 
respond to the holes not visiting by trajectories on the at- 
tractor. The north and south poles are indicated with N and 
S, respectively. The angular coordinates {9, </p) are shown for 
some point M, and axes of the Cartesian coordinates x, y, z 
are depicted 



generated by mixing of the second harmonic components 
from the oscillators on a quadrtatic nonlinear element. 
On a stage of differential rotation, we set 



d = \^a'K{^/2 - 1 - 2\/2|a|2)a, 
6= iicr7r(V2 + l-2V2|5p)6, 



(6) 



where a = ±1. In angular variables ((^, 0), these equa- 
tions reduce Xo (p = i(T7r('\/2 cos 6* -|- 1), Q — 0. Note that 
the angular velocity depends linearly on z = cos d and 
vanishes at z = — l/v^- On this stage two subsystems 
must behave like uncoupled classic non-isochronous oscil- 
lators. At small amplitudes their frequencies in respect 
to the reference point are AtJa^f, — jiaTT{\/2 ^ 1). With 
growth of the amplitudes, the oscillation frequencies un- 
dergo a shift proportional to the squared amplitude for 
both subsystems. 

Finally, on the stages of rotation an appropriate form 
of the equations is 



■jTTsb, b 



rTTSa, 



(7) 



where s = ±1. It corresponds to a conservative system of 
coupled oscillators with equal frequencies, with the cou- 
pling coefficient of such value that the energy exchange 
between the partial oscillators corresponds precisely to 
the duration of the stage. 

Now, we can write down equations for the complex am- 
plitudes embracing the complete time period T — 6. For 
this, we compose the right-hand sides as combinations of 
terms (O, dH), jT]), which are supposed to be switched 
in, or off, during the respective stages of the time evolu- 
tion. As to the terms from the equations ([T]), we account 
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them only on the stages of rotation. (Their exclusion for 
other stages is not so significant, but we do so, because it 
simplifies derivation of the Poincare map in the explicit 
form.) Finally, we arrive at the equations 



a 



2 



s^)Im{a^b*^)a- ^sb 



(\/2- 1 - 2V2|a|2)a 
■is2^(l-|a|2-|6|2)a. 



b = ie{l - 



(8) 



Im(a26*2)6 



rsa 



^icr7r(\/2 + l-2V2|6|2)6 



1^2 
2 



Here the factors a and s depend on time with period 
T = 6, and on a single period they are defined by the 
relations 



-1, l<t<2 
1, 4 < t < 5, s 
0, otherwise. 



-1, 2 < i < 3 
1, 5 < i < 6, 
0, otherwise. 



(9) 



Let us derive the Poincare map, which determines 
transformation of the state over one period T = 6 and 
describes the time evolution stroboscopically. 

Let the initial conditions for the equations (l3|) at t„ = 
nT are defined as the state vector X„ = (a„, 6„), and the 
state after a half of period is X„+i/2 = F(t,s(X„) = 

Solving equations ([8]) on each successive unit interval 
analytically, we can represent the map Fo-^^ explicitly: 



a. 



^l + (|a„|2 + |6„P)(ef-l) 



where 



a= i^7^(^/2-l-2^/2|a„|2), 



/3= ia7r(^/2 + l-2V2|6„|2), 



(10) 



D 



V2 Vl"" 



|"|b„|^-(a;f)„)" tanh(2e|a„p|fc„p) ^ ^/"^ 
|^|b„r-'-(Q„f);)^ tanh(2£|a„P|fc„P) i 



The indices cr and s become equal ±1 alternately, so the 
mapping for the complete period is defined as follows: 



— F(X„) — Fij(F_i^_i(X„)). 



(11) 



Dynamics governed by equations ([8]) or by iterations of 
the map (jlip is invariant in respect to simultaneous con- 
stant phase shift for two oscillators, i.e. to the variable 
change a ae^^ , b^be''^. Due to this, one can reduce 
the equations for two complex amplitudes to equations 
in three real variables. 

Performing the variable change ([1]), we get a set of 



differential equations 



i cr7r(zV2 + l)y - e(l - cr^ - s'^)xj/ 



y = icr7r(z\/2 + l)a;/2 + e(l - cr^ _ s'^)yx'^ 

+/is2(l - ^/3F~^]f~^Z^)y, 



(12) 



where a and s are time-dependent, as stated by formu- 
las Designating at t„ = nT the state vector as 
Xn = {xn,yn, Zji), wc cau represent the three-dimensional 
Poincare map as 



x„+i = f(x„) = fl_l(f_l,_l(x„)). 



(13) 



and the half-period map — fcr,s(x„) is expressed 

as 



Xn+l/2 



PsZn 



yn+1/2 = PQ 



-ei^l+yl)/2 sinf(z„\/2-f 1) 



+y^e«+y'^^/^ cos f(z„ 72 + 1) 



^n+l/2 — PQ 



-sx„e-"(^"+^")/2 cos f (z„V2 + 1) 



+SCTy„e'^(^"+?'")/2 sin f (z„ V2 + 1) 



where 



(14) 



P = 



{l-e-^^)^xl+yl + zl + e-^ 



^a,2e- = (=l+yS)+y2ge(x2+„2) 



The maps (fTTj) and (fT3|l are reversible. The inverse 
maps are derived from solution of the equations ([8]) or 
(I12p in the backward time; their analytic representations 
are omitted for brevity. 



III. NUMERICAL RESULTS FOR 
THE THREE-DIMENSIONAL MAP AND 
HYPERBOLIC NATURE OF THE ATTRACTOR 

In Fig. 3 a-c portraits of the attractor are shown for the 
map x„+i = f (x„) at e ~ I, fj, = 1. They are obtained by 
computation of a sufficiently large number of iterations 
after excluding the initial transient part of the orbit. As 
seen from the diagram (a), in the space (x,i/,z) the at- 
tractor is disposed on a unit sphere. In the diagram (b) 
it is represented in the angular coordinates {ip, 9), and in 
the diagram (c) as an object on the plane of the complex 
variable 



W : 



X — z + iy\/2 
x + z + ^/2 



(15) 



(a) 
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(b) 




FIG. 3: Attractor of the map (|13p at e = 1 and ^ = 1 in three-dimensional space on the unit square (a), its representation in 
the angular coordinates 6^) (b), and a portrait on the plane in stereographic projection (c). Panel (d) shows for comparison 
the portrait of the Plykin - Newhouse attractor reproduced from Ref. 24] for the Hunt model [l^]. The orientation is selected 
specially to see better the correspondence of the structure of the filaments with that in the diagram (c) 



The last corresponds to stereographic projection from 
the sphere to the plane with a use of the point C as a 
center. This point together with a neighborhood does not 
belong to the attractor of the map; so, its image occupies 
a bounded part of the plane W. 

Note specific fractal-like transverse structure of the at- 
tractor. Few initial levels of this structure are easily dis- 
tinguishable: the object looks like composed of strips, 
each of which contains narrower strips of the next level 
etc. As follows from the computations discussed below, it 
is a uniformly hyperbolic strange attractor. Its topologi- 
cal type corresponds to the Plykin - Newhouse attractor. 
The last follows from visual comparison of mutual loca- 
tion of filaments in the diagram (c) and in the Plykin 
- Newhouse attractor shown in diagram (d). (The last 
picture is taken from the paper [2J|, which reproduces 
analysis of the Hunt model [22] , definitely possessing the 
attractor of Plykin - Newhouse.) 

To compute all Lyapunov exponents for the three- 
dimensional map, joint iterations of (jl3p together with 
a collection of three equations in variations for perturba- 
tion vectors are produced. At each step, Gram-Schmidt 



process is applied to obtain an orthogonal set of vectors, 
and normalization of the vectors to a fixed constant is 
performed. Lyapunov exponents are obtained as slopes of 
the straight lines approximating the accumulating sums 
of logarithms of the norm ratios for the vectors in de- 
pendence of the number of iterations [25| . In particu- 
lar, at e = 1 and /i = 1 the Lyapunov exponents are, 
Ai = 0.9575, A2 = -1.2520, A3 = -2, and estimate of 
the attractor dimension with the Kaplan - Yorke formula 
yields = 1 + Ai/jAal « 1.765. 

To substantiate the hyperbolic nature of the attrac- 
tor, let us turn to computational verification of the 
cone criterion known from the mathematical literature 

Let us have a smooth map x = g(x) that determines 
discrete-time dynamics on an attractor A. 

The criterion requires that with appropriate selection 
of a constant 7 > 1, for any point x £ A, in the space of 
vectors of infinitesimal perturbations one can define the 
expanding and contracting cones and Cx. Here is 
a set of vectors satisfying the condition that their norms 
increase by factor 7 or more under the action of the map. 
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Cx is a set of vectors, for which the norms increase by 
factor 7 or more under the action of the inverse map 
X = g^^(x). The cones S'x and Cx must be invariant 
in the foUowing sense, (i) For any x g A the image of 
the expanding cone from the pre-image point x must be a 
subset of the expanding cone at x. (ii) For any x e A the 
pre-image of the contracting cone from the image point 
X must be a subset of the contracting cone at x. 

Let g(x) be a map corresponding to the /c-fold iteration 
of the Poincare map f (x) , where k is an integer selected 
in the course of the computations. The needed Jacobian 
matrices can be found in our case analytically, by differ- 
entiating of (fT3|) with application of the chain rule for the 
derivatives of the functional compositions. Some details 
of the computational procedure, which takes into account 
disposition of the attractor on the invariant sphere, are 
given in Appendix. 

The calculations are organized as verification of the 
required conditions for a set of point on the attractor 
obtained from multiple iterations of the map g(x). We 
check, first, the existence of nonempty expanding and 
contracting cones, and secondly, the validity of inequal- 
ities corresponding to proper inclusions of these cones. 
If these conditions are met with 7 = 1, the interval, is 
determined 7min(x) < 7 < 7max(x), in which they are 
true. 

Figure 4 shows a diagram resuming graphically results 
of verification of the cone criterion for the attractor of 
the map (jl3p at e = 1 and /i = 1. The computations 
have been performed for the map f'^(x). The diagram 
represents in logarithmic scale the values 7min(x) in gray 
and 7max(x) in black versus y coordinate of the analyzed 
points on the attractor. Observe a gray set and a black 
set, one disposed strongly above, and another strongly 
below the horizontal line 7=1. Presence of a gap sepa- 
rating these sets from the line 7 = 1 implies the positive 
result of the test. To express the result quantitatively, 
one can determine the maximum of 7min(x) and mini- 
mum of 7max(x) over the set of all processed points. As 
found, selection of the constant satisfying 0.44 < 7^ < 2.3 
ensures the required invariance of the cones. 



IV. NUMERICAL RESULTS FOR THE 
COUPLED OSCILLATORS 

In accordance with the previous section, there is a cor- 
respondence between dynamics of complex amplitudes 
in the coupled oscillators ([8]) and dynamics of the three- 
dimensional mapping (|13p possessing the hyperbolic at- 
tractor of Plykin - Newhouse. Let us illustrate with nu- 
merical results the dynamics of the coupled oscillators 
concentrating on features linked with the flow nature of 
the system, i.e. with the continuous time evolution. 

Figure 5 shows plots for the amplitudes of the coupled 
oscillators \ a\ and |6| versus time obtained from numerical 
solution of the differential equations ^ with the Runge 
- Kutta method at e = 1 and fj, = 1. Some small in 




-1 y — ► 1 



FIG. 4: A graphical illustration for verification of the hyper- 
boHc nature of the attractor for the map (|13|) . with e=l and 
fi—1. A positive result of the test follows from existence of 
the gap between the line 7=1 and the black and gray sets 
of points, which represent, respectively, the upper and lower 
edges of the intervals of 7, which met the verified conditions 



absolute value and random in phase complex amplitudes 
a and b are taken as initial conditions, so the plot depicts 
the transient process prior to the regime of chaotic self- 
oscillations. In the right-hand part of the diagram the 
dependences look like samples of a random process; that 
associates with motion on the chaotic attractor. Locally, 
some peculiarities can be seen because of the piecewise 
continuous nature of the process composed of successive 
stages. In particular, the horizontal plateaus relate to 
the stages of evolution, on which the amplitudes \a\ and 
|6| remain constant according to equations (O. Note that 
the realizations for \a\ and \b\ are interconnected: in the 
sustained regime they obey the relation |ap -I- \b\^ — 1. 

Figure 6 presents two versions of portraits of the at- 
tractor for the system ([8]) at e = 1 and n = 1. As the di- 
mension of the state space is sufficiently high (vector X = 
(a, b) is four-dimensional, and the extended state space of 
the non-autonomous system is five-dimensional), depict- 
ing the image to resolve subtle fractal transverse struc- 
ture intrinsic to the attractor is not a trivial task. For 
this, we apply presentation of the object in gray scales. 
Brighter tones correspond to pixels visiting by the rep- 
resentative point with higher probability [27| . In panel 
(a) this technique is used to draw the three-dimensional 
portrait. Angular coordinates 9) are plotted in a hori- 
zontal plane. The third variable plotted along the vertical 
axis is time, within one full period of variation of coef- 
ficients in the equations ([5]). The picture reminds rising 
and swirling smoke. In the cross-section with a horizon- 
tal plane t — (mod6) observe a fractal-like formation 
identical to the attractor of the three-dimensional map 
depicted in angular coordinates in Fig. 3 b. One more 
portrait is shown in the panel (b) on the plane of two 
values of real amplitude \a{t)\ and \a{t— 3)|, which relate 
to instants separated by a half of time period of variation 
of the coefficients in the equations ([5]). Here, again one 
can distinguish the fractal-like transverse structure linked 
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a\ 




100 200 300 t 



FIG. 5: Plots of the real amplitudes \a\ and |&| versus time in the transient process obtained from numerical solution of the 
differential equations (|8]) at e=l and 




FIG. 6: (a) Portrait of attractor for the system ((8]) in three-dimensional space {ip, 9, t). In the cross-section with the horizontal 
plane t = (mod6) observe the object identical to that shown in Fig. 3 b. (b) Portrait of the attractor in the plane of real 
amplitudes relating to one of the oscillators and separated by time interval r/2=3. Technique of representation in gray scales is 
used: brighter tones correspond to higher probability of visiting the pixels by the representative points. The parameter values 
are e=l, jj.—! 



with the dynamics of the Plykin type. This method of vi- 
sualization may be appropriate in experiments with sys- 
tems of the class under consideration. 

The Lyapunov exponents of equations ([8|) are linked 
with the exponents for the Poincare map (see (fTO| . (fTTjl ) 
by an evident relation ~ ^i/T, where T = 6 is the pe- 
riod of variation of the coefficients in the equations. The 
procedure of computation of the Lyapunov exponents Aj 
is analogous to that used for the three-dimensional map. 
Joint iterations of the map pT|) together with a collec- 



tion of four equations in variations are produced. At 
each step, Gram-Schmidt process is applied to the set of 
vectors, and normalization of them to a fixed constant 
is performed. Figure 7 present the results graphically. 
The first plot (a) shows four Lyapunov exponents Aj in 
dependence on parameter e at fixed /i = 1. In the range 
e < Ec ~ 2.03 one of the exponents is positive that means 
chaos. Among other exponents one is zero (up to nu- 
merical errors), and two are negative. Note a smooth 
dependence of the largest exponent on the parameter. 
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For larger e (strong dissipation bringing in during the 
flow down stages) chaos disappears. The second plot (b) 
shows the Lyapunov exponents versus at fixed e = 1. 
Observe that variation of /i notably effects only one expo- 
nent, which corresponds, obviously, to approach of orbits 
to the invariant sphere. Presence of a zero exponent re- 
flects invariance of the equations in respect to the overall 
phase shift. Of course, results of computations agree with 
the data from the previous section: at identical e and /i 
three nonzero exponent are equal, up to numerical errors, 
to those obtained for the three-dimensional map. 

Figure 8 shows a plot of spectral density in logarith- 
mic scale versus a frequency for a signal generated by 
one of the oscillators. It relates to regime of dynamics on 
the Plykin - Newhouse attractor interpreted in terms of 
the three-dimensional map. This spectrum is one more 
characteristic interesting in the context of possible exper- 
iments. 

In computations the standard methodic was used rec- 
ommended for non-parametric statistical estimates of the 
power spectral density. It is based on subdividing the 
whole realization on a number of parts of equal duration. 
For each part, the signal is multiplied by a smooth func- 
tion vanishing at the ends of the interval ('window'), then 
Fourier transform is applied, and finally, the squared am- 
plitudes of the frequency components are averaged over 
all the parts. A sample of the time series for the com- 
plex variable a{t) was obtained from numerical solution 
of equations ([5]) by the finite-difference method. It cor- 
responds to motion on the attractor and contains 6 • 10"' 
data points with time step At — 0.01. As seen from the 
picture, the spectrum looks continuous that corresponds 
to chaotic nature of the generated signal. The spectrum 
is almost perfectly symmetric about the reference fre- 
quency, where the spectral density is maximal. The two 
main side maximums have a level below the central one 
by about 7 dB, and their frequencies approximately cor- 
respond to the inverse value of the period of variation of 
the coefficients in equations ([S]): / w ±1/6. Apparently, 
this periodicity is a reason for the rugged form of the 
spectrum. A plot for the power spectral density for the 
second oscillator looks exactly the same because of the 
symmetry of the system. 



V. CONCLUSION AND DISCUSSION 

In the present article a system of two coupled non- 
autonomous nonlinear oscillator is introduced manifest- 
ing chaotic dynamics, which is in a direct relation with 
the concepts of the hyperbolic theory. With exclusion of 
the overall phase the Poincare map reduces to a three- 
dimensional map possessing attractor of Plykin type dis- 
posed on an invariant sphere. 

As the systems of coupled oscillators occur in many 
fields in physics and technology, it is natural to expect 
that the suggested model may be realizable. Particularly, 
it may relate to electronic devices, mechanical systems, 



objects of laser physics and nonlinear optics. The sys- 
tems with hyperbolic chaos may be of special interest for 
applications due to their robustness, or structural stabil- 
ity, that means insensitivity of the devices to variations 
of parameters, characteristics of elements, technical fluc- 
tuations, noise etc. 

Appearance of concrete examples of systems with hy- 
perbolic strange attractors makes it reasonable to apply 
for them the whole arsenal of computational methods 
accumulated in nonlinear dynamics and its applications. 
This is of evident interest both from the point of view of 
complementation of the mathematical concepts with con- 
crete and visible context (see e.g. a paper |29|), and for 
exploiting these concepts in applications. In the present 
work such results of computer studies are presented as 
realizations, attractor portraits, Lyapunov exponents, es- 
timate of dimension, power spectral density. 

It is worth mentioning some possible modiflcations of 
the model. 

• It is easy to suggest a version of the equations, 
in which the temporal dependence of the coeffi- 
cients will be piecewise smooth. For this, one can 
introduce in the equations x = f (x, i) a smooth 
common time-dependent factor vanishing at the 
junctions of the stages, integral of which over a 
stage duration equals 1. An appropriate variant is 
X — (2 sin^ 7rt)f (x, i). The Poincare map remains 
the same, and the nature of the attractor is not 
changed too. (Hunt used a similar trick in his the- 
sis [2a|.) 

• As mentioned, the "self-oscillatory" term in the 
equations proportional to /i may be retained on all 
stages of the dynamics. 

• Working with a version of the model describing by 
the three-dimensional set of equations, it is possible 
to simplify the form of the nonlinearity excluding 
the operation of extracting the square root, and 
setting the respective term to be ^lil — x^ — y^ — z^). 
This modification does not influence dynamics on 
the attractor belonging to the invariant set + 

= 1. (In amplitude equations for a and b such 
modification leads rather to complication because 
of increase of degree of the nonlinear terms.) 

• Taking into account structural stability of the hy- 
perbolic attractor, many other modifications can 
be done, which do not change the nature of the 
attractor, while the variations are not too large. 
In particular, it is possible to introduce a model 
with smooth analytical variation of the coefficients 
in time in the non-autonomous equations psj . 

Formally, in our complex amplitude equations the at- 
tractor should not be interpreted as uniformly hyper- 
bolic, because of presence of a neutral direction in the 
state space associated with the overall phase. Instead, 
it has to be related to the class of partially hyperbolic 
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FIG. 7: Four Lyapunov exponents versus parameter e at /i = 1 (a) and versus parameter /i at e = 1 (b) for the Poincare map 
represented in terms of complex amplitudes. A zero exponent occurs due to invariance of the equations in respect to the overall 
phase shift 




FIG. 8; The power spectral density for one of the coupled 
oscillators versus the frequency computed by processing a re- 
alization obtained from numerical integration of the equations 
(131) at e = 1 and n — 1 



attractors [a, l23|. Nevertheless, in the form used here in- 
variance of the equations in respect to the overall phase 
is exact; it means that one can accept a rightful agree- 
ment not to distinguish states distinct only in the phase, 
and, in this sense, treat the dynamics as true hyperbolic. 

However, in systems, for which description in terms of 
slow complex amplitudes will be an approximation, one 
can expect appearance of peculiarities associated with 
features of the partially hyperbolic attractor. If the 
deflections from the slow-amplitude approximation are 
small, the overall phase will manifest slow random walk, 
while the dynamics of the rest variables will retain its 
character because of intrinsic robustness. 

In dependence on parameters, the suggested model can 
manifest chaos and regular (periodic) dynamics. So, it 
may serve as an object for principal and interesting stud- 
ies of scenarios of the onset of hyperbolic strange at- 
tractors in the course of parameter variations (see e.g. 
[sol [3l| ) . Insufficient progress in this research direction 
may be explained particularly by the fact that no realistic 
examples of concrete systems undergoing such transitions 
were known. 

The work was performed under support of RFBR - 
DFG grant No 08-02-91963. 



Appendix 

For a three-dimensional dissipative map x = g(x), 
X, X e let us consider the procedure for verification 
of the cone criterion, bearing in mind the situation that 
one expanding and two contracting directions present in 
the state space, and the attractor is located on the in- 
variant sphere. The map is supposed to be reversible: 
any state vector x has a unique image x = g(x) and a 
unique pre-image x = g^^(x). 

Let the derivative matrix of the map g at x be v = 
dx(g), which acts in the tangent space of vectors u = 
{ui,U2,u^}- Via the auxiliary symmetric matrix b = 
v'^v, where the superscript T means the transpose, the 
norm of the vector u = vu is expressed as 



u^bu. 



(16) 



The expanding cone at the point x is a set of vectors 



5x = 



|u u^bu>7^u^u| 



(17) 



With the same matrix b one can define a pre-image of the 
contracting cone relating to the point x — g(x), namely. 



C4 = |u 72u'^bu < u'^uj 



(18) 



Now, let us consider an inverse map x = g^^(x) and 
its matrix derivative w = dx(g^^). Via the auxiliary 
symmetric matrix a = w'^w we represent the norm of 
the vector u = wu as 



- I|2 
U 



T - 

u au. 



(19) 



With the matrix a we define the contracting cone at x 
as a set 



Cx = {u I u^au > 



(20) 
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and a cone that is an image of the expanding cone relating 
to the point x — g~^(x), namely, 



54 = {u I -f'^u^au < u^u } 



(21) 



In computations it is necessary to check, first, existence 
of nonempty cones satisfying the definitions, and, second, 
validity of the inclusions S"- C and C4 C Cx- 

As the attractor is placed on an invariant sphere , let 
us assume that x S S*^ . To define a convenient orthogonal 
basis {ii, '12, is} we take as is a unit vector directed along 
the radius at x, and the unit vectors ii, 12 are taken in 
the tangent plane. To be concrete, we require the matrix 
element 612 to vanish, and the inequality bn > 622 to 
hold. 

The conditions of required inclusion of the cones are 
formulated in terms of quadratic forms associated with 
the matrices b = b — F^e and a = a — F^e, where e is 
the unite matrix. A constant factor F is assumed to be 
equal 7, or I/7, considering the expanding, or contracting 
cones, respectively. Note that the matrices a and b are 
symmetric: bij — bji, Uij = Uji. 

Equations 

biiul + b22ul + &33M3 + 2613U1M3 + 2&23W2U3 = (22) 
and 

aiiul+a22ul+a33ul+2ai2UiU2+2ai3UiU3+2a23U2U3 = 

(23) 

determine the borders of the cones. By variable change 
u[ = ui + 6ii^6i3W3, "2 = "2 + 622^623^3, "3 = ""3 (24) 



the quadratic form in 



is reduced to a standard form: 

(25) 



^'iiu'i + b22u'2 + b'^^u'l — 



while the equation ([23]) becomes 

aiiu I + a22U 2 + 0,33^ 3+ 
+2ai2u' lu' 2 + 2ai3u'iu'3 + 2a'23u'2u' 3 = 0. 



Here 



and 



b' 



33 



^33 



6n&?3 



"22 "23 



(26) 



(27) 



-^33 



033 



: ai3 - aiib^lbi3 - 012622^623, 

= 023 - 0126^1^613 - 022622^623, 

aii6n^6f3 + 022622^6^3 - 2ai3b^^bi3 
-2023622^623 + 201261^1^613622^623. 



'13 



*23 



(28) 



In the cross-section by a plane u'l — const the equa- 
tions (|25|1 and ([26|) determine some curves of the second 
order; their types and mutual location have to be re- 
vealed in the course of computations. To have situation 



of inclusion required by the criterion, these curves must 
be ellipses. 

First, in computations we check the inequalities 611 > 
0, 622 < 0, 633 < 0. If they are true, the equation (^5]) 
defines an ellipse. 

To determine the type of the curve given by (I26p , we 
compute the invariants 



/ = 022 + O33, D = 



022 O23 

023 O33 



011 O12 O13 

012 022 O23 

013 O23 O33 



(29) 

and check that I? > 0, and A/I > {]. Then, in accordance 
with the theory of conic sections, the equation ((26|) also 
defines an ellipse. 

Let us formulate a convenient and simple sufficient con- 
dition of location of the second ellipse inside the first one. 
Renormalizing variables 



i = yb^iu'i, 77 = V-^22W2, C = \/-~hW3 (30) 

and setting u'l = 1/ \/6ii, transforms the first ellipse to a 
unit circle 



if+e = i. 



Then, the equation for the second ellipse is 



_£22 2 
022 ' 



2ai2 



2°13 A I ail — n 



Qll _ I 



(31) 



(32) 



and its center is located at 

•no = 

Ql3Q22-t-Ql2Q: 



0220^3-0^3^ 



622. 
fell ' 



(33) 



Co 



a22 0>-:i'i—0.2 



"bll 



In variables 77 = rj—rjo, ( = C^Co the equation becomes 



022 -2 
622 



2a'n 



23 



V0226 



(34) 



33 



^^33 



where 



i?2 



Oil I 022 J.2 
bii ^ 622 '0 



V-biib22^'^ ^-biib'^^ 

Computing the lesser root of the square equation 

l'^ + (022/622 + 033/6^53)? -I- (022033 - 4a'23^)/(622633) = 0, 

(36) 

one finds out the semi-major axis R/ ^/Imin- A sufficient 
condition for the ellipse to be located inside the unit disc 
is inequality 



Y^>22bi 

2o', 



(35) 



< 1. 



(37) 
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If all the named conditions are true at F = 7 and at 
r = 1/7, one can deduce about the correct inclusion 
for the expanding and contracting cones at the analyzed 
point X. 

Indeed, in variables 77, C the cross-section of the ex- 
panding cone iSx is the closed unit disc, and cross-section 
of the cone is represented by the closure of interior 
of the small ellipse obtained at F = 7, so, the required 
inclusion C is valid (Fig. 9 a). On the other hand, 
cross-section of the contracting cone Cx is a closure of 
exterior of the small ellipse obtained at F = I/7. Cross- 
section of the cone corresponds to a closure of exterior 
of the unit circle. Hence, the inclusion C Cx is valid 
(Fig. 9 b). 



r=Y r=i/Y 




(a) (b) 

FIG. 9; Cross-sections of the cones for the case of three- 
dimensional map with one expanding and two contracting 
directions. A picture of proper inclusion is shown (a) for ex- 
panding cones C S'x and (b) for contracting cones C4 C Cx 
A circle circumscribed around the ellipse is located inside the 
unit disc if the condition (19) is valid 



Computations in the present work were organized as 
verification of the cone criterion for a set of points on the 
attractor obtained from multiple iterations of the map 
g(x). The inclusions were checked with a help of the 
inequality p7|) . 

Because of smoothness of the map under study, the ob- 
jects considered in the context of the cone criterion (ma- 
trices, quadratic forms and their invariants) depend on 
the state variables in smooth manner, as they are deter- 
mined by dynamics on finite time intervals. As follows, 
validity of the conditions at some point x with distant 
from 1 constant 7 implies that the cone criterion holds as 
well in a neighborhood of x (as wider, as larger the value 
I7 — 1| is). A positive result of the test for a represen- 
tative set of points implies validity of the cone criterion 
on the whole attractor, if it is completely covered by the 
union of the mentioned neighborhoods. Practically, such 
situation is achieved by increase of the number of itera- 
tions and, respectively, a number of tested points on the 
attractor. 

It is convenient not to fix in advance the constant 7, 
but to arrange the computations as follows. First, at 
each point x we compute the matrices a and b and check 
all the formulated conditions for 7=1. If they hold, the 
program determines an allowable interval of 7. For this, 
the program simply enumerates the 7 values with a small 
step in a sufficiently wide range. Attractor is recognized 
as hyperbolic if a gap of a finite width separates the ob- 
tained sets of top and bottom edges of the intervals from 
the axis 7 = 1 on the plot of 7 versus some dynamical 
variable characterizing location of the analyzed point x. 
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